module mWigner

        use mPrecision

        implicit none

        interface wigner
                module procedure wigner_integer
                module procedure wigner_D
        end interface wigner

        private
        public :: wigner

contains
        
        real(D) function wigner_integer(l1, l2, l3)

                use mMisc

                integer, intent(in) :: l1, l2, l3

                interface
                        function gamma(x)
                                real(8) :: gamma
                                !DEC$ ATTRIBUTES C,ALIAS:'gamma' :: gamma
                                real(8), intent(in) :: x
                        end function gamma
                end interface
#ifdef PGF90
!$PRAGMA C(gamma)
#endif ! PGF90

                integer :: l, lo2, lo2p1, lp1

                l = l1 + l2 + l3


                if( odd(l) ) then
                        wigner_integer = 0.0_D
                else
                        lo2 = l / 2
                        lo2p1 = lo2 + 1
                        lp1 = l + 1
                        wigner_integer = real(gamma( real(lo2p1, 8) ) &
                                            - gamma( real(lo2p1 - l1, 8) ) &
                                            - gamma( real(lo2p1 - l2, 8) ) &
                                            - gamma( real(lo2p1 - l3, 8) ) &
                                       + 0.5_8 * (gamma( real(lp1 - 2 * l1, 8) ) &
                                            + gamma( real(lp1 - 2 * l2, 8) ) &
                                            + gamma( real(lp1 - 2 * l3, 8) ) &
                                            - gamma( real(l + 2, 8) )), D)
                        wigner_integer = exp(wigner_integer)
                        wigner_integer = wigner_integer * (-1) ** lo2
                end if

        end function wigner_integer

        real(D) function wigner_D(l1, l2, l3)

                use mMisc

                real(D), intent(in) :: l1, l2, l3

                real(D) :: l, lo2, lo2p1, lp1

                interface
                        function gamma(x)
                                real(8) :: gamma
                                !DEC$ ATTRIBUTES C,ALIAS:'gamma' :: gamma
                                real(8), intent(in) :: x
                        end function gamma
                end interface
#ifdef PGF90
!$PRAGMA C(gamma)
#endif ! PGF90

                l = l1 + l2 + l3

                lo2 = l / 2.0_D
                lo2p1 = lo2 + 1.0_D
                lp1 = l + 1.0_D
                wigner_D = real(gamma( real(lo2p1, 8) ) &
                              - gamma( real(lo2p1 - l1, 8) ) &
                              - gamma( real(lo2p1 - l2, 8) ) &
                              - gamma( real(lo2p1 - l3, 8) ) &
                         + 0.5_8 * (gamma( real(lp1 - 2 * l1, 8) ) &
                                  + gamma( real(lp1 - 2 * l2, 8) ) &
                                  + gamma( real(lp1 - 2 * l3, 8) ) &
                                  - gamma( real(l + 2, 8) )), &
                               D)
                wigner_D = exp(wigner_D)

        end function wigner_D

end module mWigner
